misclass_phi1_kpred=function (form,  bmat = 0, print.summary = TRUE, 
                        data = NULL) # Here form must be y_both~x1+x2+x3+...xp+z1+z2....
{
  probit <- glm(form, family = binomial(link = "probit"), data = data) 
  #cat("Standard Probit Estimates", "\n")
  #print(summary(probit))
  ysum <- model.frame(form, data = data)[, 1]
  xzmat <- model.matrix(form, data = data) # Design matrix (1,x1+x2+x3+...xp,z1,z2)
  nxz = ncol(xzmat)  
  probitX <- glm(ysum~ xzmat[,2:(nxz-2)], family = binomial(link = "probit"), data = data)  # fit probit  considering just  x1+x2+x3+...xp
  xmat<- model.matrix(ysum~ xzmat[,2:(nxz-2)]) # Design matrix (1,x1+x2+x3+...xp)
  if (identical(bmat,0)){
    nbxzmat <- length(probit$coef)
    bxmat <- probitX$coef
    bxzmat <- probit$coef
    bmat<-c(bxmat,bxzmat)
  }
  nx = ncol(xmat)
  
  nk = length(probit$coef)+ length(probitX$coef)
  bstart <-  bmat
  # names(bstart)<-c("beta0","beta1",....,"betap","eta0","eta1","eta2","eta3")
  nameb=rep("beta",nx)
  namee=rep("eta",nxz)
  namebeta=outer(nameb, seq(0,(nx-1)), function(x,y) paste(x,y,sep=""))[1,]
  nameeta=outer(namee, seq(0,(nxz-1)), function(x,y) paste(x,y,sep=""))[1,]
  namebstar<-c(namebeta,nameeta) 
  names(bstart)<-namebstar
  
  logl <- function(b) {
    eysum<-(1-pnorm(xzmat%*% b[(nx+1):nk]))*pnorm(xmat%*%b[1:nx]) 
    out<--sum(ifelse(ysum==1,log(eysum),log(1-eysum)))
    return(out)
  }
  
  nlfit<-optim(bstart,logl, gr=NULL, method="BFGS",control=list(maxit=500), hessian=TRUE) # the function probit.lf is defined above
  # For default maxit=100 and the iteration limit may need to be larger to get convergence. 
  vmat <- solve(nlfit$hessian)
  stderr <- sqrt(abs(diag(vmat)))
  tval <- nlfit$par/stderr
  pval <- 2 * (1 - pnorm(abs(tval)))
  
  cname <- c("Estimate", "Std. Error", "z-value", "Pr(>|z|)")
  if (print.summary == TRUE) {
    outmat <- cbind(nlfit$par, stderr, tval, pval)
    #rownames(outmat) <- names(probit$coef)
    colnames(outmat) <- cname
    print(round(outmat, 5))
    cat(" ", "\n")
  }
  
  out1 <- list(nlfit$par, stderr, vmat, nlfit$convergence, 
               -nlfit$value)
  names(out1) <- c( "estimate", "stderr", "vmat","convergence", "minimum")
  return(out1)  
}
